Thermal phase transitions in Artificial Spin-Ice 
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We use the sixteen vertex model to describe bi-dimensional artificial spin ice (ASI). We find excellent agree- 
ment between vertex densities in fifteen differently grown samples and the predictions of the model. Our results 
demonstrate that the samples are in usual thermal equilibrium away from a critical point separating a disordered 
and an anti-ferromagnetic phase in the model. The second-order phase transition that we predict suggests that 
the spatial arrangement of vertices in near-critical ASI should be studied in more detail in order to verify whether 
they show the expected space and time long-range correlations. 
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Hard local constraints produce a rich variety of collec- 
tive behavior such as the splitting of phase space into dif- 
ferent topological sectors and the existence of "topological 
phases" that cannot be described with conventional order pa- 
rameters Q. In geometrically constrained magnets, the lo- 
cal minimization of the interaction energy on a frustrated unit 
gives rise to a macroscopic degeneracy of the ground state 0, 
unconventional phase transitions EH), the emergence of a 
"Coulomb" phase with long-range correlations 0[6l and slow 
dynamics OH) in both 2D and 3D systems. In this work we 
focus on a paradigm with these features: spin-ice, a class of 
magnets frustrated by the ice-rules I9UTT1. 

In natural spin-ice the ice rules are due to two facts: the 
pyrochlore lattice structure, in which rare earth magnetic ions 
sit on the vertices of corner sharing tetrahedra, and the ferro- 
magnetic interaction between the Ising-like moments that are 
forced to lie on the edges joining the centers of neighboring 
tetrahedra by crystal fields. The energy of each unit cell is 
thus minimized by configurations with two spins pointing in 
and two out of the center of the cell. All configurations sat- 
isfying the local constraint are degenerate ground states if the 
ice-rule preserving vertices are equally probable. This is the 
same mechanism whereby water ice has a non-vanishing zero- 
point entropy that is, indeed, remarkably close to the one of 
the Dy2Ti2C>7 spin-ice compound lfT2l . These materials re- 
ceived a renewed interest in recent years when a formal map- 
ping to magnetostatics suggested to interpret the local config- 
urations violating the spin-ice rule as magnetic charges lfT3l . 
However, the detailed study of such defects in 3D remains a 
hard experimental task. 

Bi-dimensional Ising-like ice-models had no experimental 
counterpart until recently when it became possible to manu- 
facture artificial samples made of arrays of elongated single- 
domain ferromagnetic nano-islands frustrated by dipolar in- 
teractions. The beauty of artificial spin-ice (ASI) is that the 
interaction parameters can be precisely engineered — by tun- 
ing the distance between islands, i.e. the lattice constant, or by 
applying external fields — and the microstates can be directly 
imaged with magnetic force microscopy (MFM) lfT4l . These 
systems should set into different phases depending on the ex- 
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Figure 1 . Sixteen vertex configurations and their statistical weights 
uji oc exp (—(3 a). The first six- vertices, vi, . . . ve, verify the ice- 
rule and have vanishing magnetic charge. Vertices v± , . . . , V4 carry 
a dipolar moment while v$ and vq do not. The last ten vertices are 
"defects": vj and vg have magnetic charge ±4 and no dipole mo- 
ment while vq, . . . , vie have magnetic charge ±2 and a net dipole 
moment. 

perimental conditions fT5ll . However, the lack of thermal fluc- 
tuations due to the high energy barriers for single spin-flips, 
had prevented the observation of the expected two-fold de- 
generate antiferromagnetic (AF) ground state. Recently, these 
problems have been partially overcome via (i) the gradual 
magneto-fluidization of an initially polarized state lfT6l and 
(ii) the thermalization of the system during the slow growth of 
the samples ifTTl . Although with de-magnetization the actual 
AF state of square ASI was never reached, the statistical study 
of a large number of frozen configurations of samples with up 
to 10 6 vertices at interaction-dominated low energies became 
possible lfT6l . With the procedure proposed in (ii) large re- 
gions with AF order were formed in a few samples when suf- 
ficiently small lattice constant and weak disorder were used. 
Whether the sampling of a conventional thermal equilibrium 
ensemble is achieved in this way is a question that was raised 
in (171 [18) and that we will address here. 

The purpose of this letter is to interpret and explain very re- 
cent experimental observations on ASI 1 16 - 18 ]. With this aim 
we consider a simple schematic model for 2D dipolar spin- 
ice, the sixte en-vertex model (see Fig. [I]) [4, 8], where dipolar 
interactions beyond nearest-neighbor vertices are neglected. 
We compare the equilibrium and out-of-equilibrium proper- 
ties of the model with the experimental results and we match 
the behavior of several observables (such as the densities of 
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Figure 2. (Color online.) Characteristic configurations of the ground 
state (d = e = 0) of three different phases of the sixteen-vertex 
model on the square lattice, (a) Antiferromagnetic (c-AF) order. 
Reversing the central (red) loop yields an elementary excitation, 
(b) Collective spin liquid (SL) phase, (c) Frozen ferromagnetic (FM) 
order. The reversed vertical string (red line) is an extended excita- 
tion. 

different vertex types (n^)) with data on ASI samples fT7l[T8ll . 
This approach is at face value similar but actually very differ- 
ent from the one used in previous studies llUOS), wriere data 
were fitted in terms of single independent vertices (with no 
topological interactions), as (n^) ex exp(— /3 e ff e^), with the ef- 
fective temperature /3 e ff introduced as a fitting parameter. On 
the contrary, in our approach frustrated interactions are fully 
taken into account. This analysis allows us to address the issue 
of thermalization of ASI samples and to make several predic- 
tions that could be tested in the lab. 

The samples and the model In their simplest setting ASI 
are 2D arrays of elongated single-domain permalloy islands 
whose shape anisotropy defines Ising-like spins arranged 
along the edges of a regular square lattice. Spins interact 
through dipolar exchanges and the dominant contributions are 
the ones between neighboring islands across a given vertex. 
The sample is frustrated since no configuration of the sur- 
rounding spins can minimize all pair- wise dipole-dipole inter- 
actions on a vertex. In samples with no height offset (h = 0) 
2D square symmetry defines four relevant vertex types of in- 
creasing energy, where the c vertices (see Fig. [I]) take the low- 
est value, leading to a ground state with staggered c-AF order 
[see Fig.[2ja)]. Note that the relative energies of the different 
vertex configurations could be tuned by h in such a way that 
the ground state displayed FM order [see Fig.[2|c)] fl5l . 

We mimic the experimental samples with the sixteen-vertex 
model defined as follows: Ising spins sit along the edges 
of an L x L square lattice. Long-range interactions beyond 
next nearest neighbor spins are neglected and the energies of 
the sixteen vertex configurations are attributed as explained 
above. The total energy is given in terms of populations, n^, 
of distinct vertex types, E = Y^\=i n % e %-> eacn w i m a given 
magnetostatic energy . Note that, despite the simple form of 
E, nontrivial frustrated topological interactions between ver- 
tices are present, due to the fact that each pair of neighboring 
vertices shares the spin sitting along the edge joining them. 
The model is coupled to a thermal bath at inverse tempera- 
ture p. We introduce the statistical weights uoi ex exp(— /3e^), 
that in the literature of vertex models are usually referred as 

a = UJi = Co>2, b = UJs = Co>4, C — UJr> — Co>6, d = Co>7 = Co>8 

and e = ujg = . . . = ujiq fU (see Fig. [TJ. The equilibrium 
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Figure 3. (Color online.) Average densities of different vertex types 
as a function of /3eo. Full symbols with error bars are experimental 
data tl8l . (ni) e x P - Empty symbols with dotted lines correspond to 
the equilibrium CTMC data, (ni) S im- The CVBP analytic solution, 
(ni)MF, of the sixteen-vertex model is shown in solid lines. 

properties of the model can thus be described in terms of four 
different statistical weights satisfying c > (a — b) > e > d. 

We study the equilibrium and out-of-equilibrium properties 
of the sixteen-vertex model in two ways. We perform numer- 
ical simulations of the 2D model with a Monte-Carlo algo- 
rithm with single- spin updates improved with the rejection- 
free Continuous Time set-up (CTMC). We also employ an 
analytic approach put forward in fl9l , based on a sophisti- 
cated version of a cluster variational mean-field Bethe-Peierls 
(CVBP) formalism, by defining the model on a coordination- 
four tree made of square plaquettes (an extension over the 
tree of single vertices is necessary to correctly describe an AF 
phase populated by finite loop fluctuations, see Fig. [2]). De- 
tails of the calculations have already been extensively reported 
in fT9lL where the phase diagram was derived for generic val- 
ues of the ratios a/c, b/c, d/c, e/c, showing that the CVBP 
technique describes with extremely good accuracy the equi- 
librium properties as compared to MC simulations. 

Density of vertices. In the experiments in [Q31 El the 
thickness of the magnetic islands grows by deposition (at con- 
stant temperature and all other external parameters within ex- 
perimental accuracy) on fifteen lattices with five different lat- 
tice constants and using three material underlayers (Si, Ti, 
Cr) in each of them. (We refer to the supplementary material 
in llT8l for more details on the parameters and materials used.) 
The Ising spins flip by thermal fluctuations during the growth 
process. However, as the energy barrier for single spin flips 
increases with the size of the islands, once a certain thickness 
is reached the barrier crosses over the thermal fluctuations of 
the bath, /c#T, and the spins freeze. At the end of the growth 
process, the spin configurations are imaged with MFM and 
the number of vertices of each kind are counted on five in- 
dependent square areas with, roughly, 27 x 27 vertices each. 
Average values (and statistical errors) are also estimated. 

Our analysis is as follows. For each set of vertex concen- 
trations measured in f\M at a given lattice constant, we de- 
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Figure 4. (Color online.) Staggered magnetization M_ (a) and spe- 
cific heat C (b) as a function of the distance to the inverse critical 
temperature (3 C for system sizes L = 30, 40, 50, 60. The (red) 
solid lines are the results of the CVBP analytic calculation. 

termine the statistical weights c, a, e, and d — and thus f3e c , 
f3e a , f3e e , and f3cd — that better match the data, by solving 
the sixteen- vertex model with the CVBP approximation and 
imposing: 



(fH)e 



( n i) 



MF 



J2 C me- n ^ d ^ 



(i) 



where the sum runs over all possible vertex configurations C. 
Surprisingly enough, it turns out that the energy ratios are ap- 
proximatively constant within the statistical errors for almost 
all the lattice spacings of the as-grown samples, and coincide 
with the values used by Nisoli et ah fl6l . Such energy ra- 
tios can be rationalized in term of magnetostatic exchanges 
due to dipolar interactions between the islands of the sample 
associated to each vertex configuration. More precisely, each 
dipole is considered as a pair of oppositely charged monopoles 
sitting on the vertices, and only Coulomb interactions be- 
tween monopoles around a single vertex are taken into ac- 
count, yielding: e c = 2(1 - 2y/2)/£, e a = -2/1, e e = 0, 
€d = 2(2y/2 + l)/i(i being the lattice constant). As a conse- 
quence, the statistical weights of the sixteen-vertex model can 
be expressed in term of a single energy scale eo as c = 1, a = 
e _/3e °, e = e -/3re€ °, and d = e _/3rde °, with the energy ratios 
r e and r d being equal to r e = (2y/2 - l)/(2y/2 - 2) ~ 2.207 
and r d = 2\/2/(\/2 - 1) ~ 6.828. Remarkably, no fitting 
parameter nor effective temperature needs to be introduced to 
describe the data: f3 is the true inverse temperature at which 
experiments are performed and eo is an energy scale that only 
depends on the lattice constant and the underlayer, and which 
could be in principle determined from a microscopic calcula- 
tion. 

In Fig. [3] we plot the thermal averages of the vertex densi- 
ties as a function of the energy scale fie^. The CTMC results 
((ni)sim> empty symbols) and CVBP calculations ((ti^mf, 
solid lines) are in excellent agreement. They show a second 
order phase transition from a paramagnetic (PM) phase at low 
P to an AF phase dominated by type-c vertices at high /3, 
where both translational invariance and spin reversal symme- 
try are spontaneously broken. We also include data from |[T8l 
(full symbols) that are in remarkable good agreement with the 
sixteen vertex model curves. (Except for a few points that turn 
out to be quite close to the critical temperature.) 

Such a good agreement confirms that dipolar interactions 
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Figure 5. (Color online.) Space-time correlations C(r, t) after a 
quench from ft = to fi c « 1.2 (a) and (3 — 1.36 (b), in a system 
with L = 60. The colored lines-points are data taken at different 
times. The dotted black lines are the t — » oc equilibrium functions, 

C(r,t -> oo) ~ r' 2 ^' 71 with 77 = 1/4, and C(r,t -> 00) ~ 



A exp(-r/£) + Ml, with Ml w 0.73, f « 3. The inset shows the 
data in panel (b) where the r-axis has been rescaled by t 
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beyond nearest neighbor vertices do not play a prominent role 
and that the sixteen-vertex model mimics ASI samples. It also 
points out that the simple picture put forward in [16 ] provides 
a good estimate of the energy ratios between different vertex 
types. Finally, and most importantly, these results strongly 
suggest that the gradual growth of magnetic islands (171 [HI 
seems to sample the conventional thermal equilibrium ensem- 
ble for most of the experimental points. 

Phases and phase transition. In the following we focus 
on the PM/AF phase transition. In Fig. [4] we present equi- 
librium CTMC data for: (a) the order parameter describing 
c-AF ordering, M_ = \({\m x _\) + (\m y _\)) where m x _l y are 
the staggered magnetizations along the horizontal and vertical 
directions; (b) the heat capacity C = L~ 2 ((E 2 ) — (E) 2 ) as 
a function of the distance to the critical inverse temperature, 
(/? — P c ) I ' j3 c . The data clearly show the presence of a second- 
order phase transition from a conventional PM phase to a stag- 
gered AF phase as fj is increased above f3 c eo. The panels dis- 
play the analytic results with solid red lines that yield a sys- 
tematic shift of the critical point by about 10% towards higher 
temperature, as expected for mean-field calculations. We find 
f3^ F = 1.05 and (3 s c im = 1.2, where here and in the follow- 
ing we measure f3 in units of eo . 

We determine j3 c = 1.204 ± 0.008 independently with a 
non-equilibrium relaxation method l20l , by identifying f3 c as 
the inverse temperature at which the staggered magnetization 
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Figure 6. (Color online.) Snapshots of L = 10 2 samples at t c± 
10 2 , 10 3 , 10 4 MCs after a quench from /5 = to the critical point, 
/3 C = 1.2 (a) and into the AF phase, /3 = 1.36 (b). Pink regions are 
c-AF ordered, green points correspond to a, b-FM vertices, red and 
blue points correspond to oppositely charged defects of type e. 

has an algebraic decay, M_(t) ~ ^-^/O^) as a function of 
time, where fj and ^ are the equilibrium critical exponents as- 
sociated to the order parameter and the correlation length, re- 
spectively, and z c is the dynamical exponent. Away from this 
temperature an exponential decay is instead observed. This 
confirms that the criticality of the SL phase (U [2T1 is broken 
by the presence of defects at finite temperature. 

Correlation functions. Stochastic thermal evolution and 
canonical equilibration not only yield the vertex densi- 
ties but also the correlation functions. Figure [5ja) shows 
the space-dependence of the two-point function C(r, £) = 
L~ 2 (^2 i • Sij(i)Si+r,j+r(t)) at different times after a 
quench from /3 = to /3 C = 1.2 (log-log scale). Sij denote 
the spins sitting along the edges of the L x L square lattice, 
with Si j — +1 if the spin points right or up and Sij — — 1 
otherwise, r is measured in units of t/y/2. These data have 
been averaged over 10 3 independent runs of a system with 
L = 60. At large times the curves approach the equilibrium 
asymptotic law characterized by C(r,t —> oo) ~ r~ v at the 
critical point, with an exponential cut-off. As shown in the fig- 
ure (black dotted line), the exponent r] remains equal to 1/4 in 
the sixteen-vertex model as it was argued in fT9) . This value 
coincides with the exact result of the eight- vertex (H ED and 
2D Ising models. Note that the numerical data show clear 
signs of finite size effects when r w L. 

In Fig.^b) we show the behavior of C(r, t) after a quench 
into the c-AF phase (/3 = 1.36). As shown in E3, the 
approach to equilibrium is fast if the initial state is a T = 
ground state whereas it is very slow and occurs via a coars- 
ening process if the initial condition is a disordered high tem- 
perature one, as for the curves shown in the figure. In equi- 
librium correlations decay exponentially as C(r,t — » oo) ~ 
A exp(— r I £) + M 2 (dotted black line). The asymptotic value 
M 2 w 0.73 is consistent with the equilibrium staggered 



magnetization shown in Fig. [4j and the correlation length is 
£(/? = 1.36) « 3. On the other hand, out-of-equilibrium spa- 
tial correlations decay to zero at large distances. The data for 
C (r, t) at different times shown in Fig. |5jb) collapse onto a 
single curve when the length variable r is rescaled by t 1 ^ 2 
(see the inset). This scaling is accurate in a certain time 
window (i.e. the coarsening regime): it fails at times larger 
than « 2600 MCs (as shown in the inset) and smaller than 
w 100 MCs. The snapshots in Fig. [6] and the t 1 ! 2 scaling 
strongly suggests that the system follows a curvature driven 
type of dynamics during the coarsening regime |23) . 

These results imply that the samples obtained by using 
the rotating field protocol fT6lL where no correlations beyond 
first-neighbors were observed, have not achieved equilibrium. 
On the contrary, the as-grown samples in (T71 [18) are likely 
to be near equilibrium. However, as also shown by our nu- 
merical results, close to the phase transition critical slowing 
down sets in, and in the whole AF phase slow coarsening dy- 
namics emerge. It remains therefore to be understood whether 
critical and subcritical samples have/w/fy equilibrated. To set- 
tle this issue one should grow samples with a slower deposi- 
tion rate and measure, if possible, time-dependent observables 
during growth (such as the staggered magnetization and two- 
times correlation functions). Another possibility is to analyze 
the spatial correlations. Indeed, the facility of imaging the 
microstates both numerically and experimentally makes such 
study very appealing. As an example, in Fig. [6] we show snap- 
shots of microscopic configurations of L = 10 2 systems after 
quenches from j3 = to j3 c = 1.2 (a) and f3 = 1.36 (b). The 
first two panels in each row are out-of-equilibrium while the 
last ones show typical equilibrium configurations at the crit- 
ical point and inside the AF phase respectively, cfr. Fig. [5] 
Large domains of ground state AF order form, separated by 
domain walls made by a and b vertices. The size of the AF 
domains increases with time as equilibrium is approached. 
These snapshots could be compared with MFM images of ASI 
configurations of H7][T8) at the corresponding values of (j (the 
experimental points closest to f3 c = 1.2 and (j = 1.36 cor- 
respond to the samples produced using a Ti underlayer with 
lattice constants i = 466 nm and t = 433 nm, respectively). 

Conclusion. In this letter we studied the equilibrium and 
out-of-equilibrium properties of the 2D sixteen-vertex model 
using MC simulations and a sophisticate cluster variational 
Bethe-Peierls approach, and we compared our results to the 
recent experimental data on ASI fT6lfT8) . We showed that the 
model describes with very good accuracy the behavior of the 
densities of the different vertex types of the ASI samples ob- 
tained by gradual deposition of magnetic material on a square 
pattern as the lattice constant and the underlayer disorder are 
changed, resulting in a change of the statistical weights of 
the different vertex types. This implies that the experimental 
samples of [031 [18) are at — or at least very close to — thermal 
equilibrium. It is important to point out that our interpretation 
does not require any fitting parameter as the effective tem- 
perature introduced in (l6l[T8). We reveal the presence of a 
second order phase transition from a conventional high tern- 
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perature (large lattice constant, strong disorder) PM phase to a 
low temperature (small lattice constant, weak disorder) stag- 
gered AF phase. Such a phase transition could have a major 
impact on thermalization and full equilibration, and could be 
unveiled by measuring long-range spatial correlations in the 
experiments. 

We close by insisting upon the fact that, although vertex 
models avoid all the complications of (long-range) dipolar 
interactions, they provide a very good schematic framework 
to study ASI from a theoretic perspective. The excitement 
around these samples as well as the intriguing excitation prop- 
erties of spin-ice (emergence of magnetic monopoles and at- 
tached Dirac strings fT3ll ) should encourage their study from 
a novel and more phenomenological perspective. 
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